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Abstract 

We are concerned with central differencing schemes for solving scalar hyperbolic con- 



servation laws. We compare the Kurganov-Tadmor (KT) two-dimensional (Kurganov 
and Tadmor , 2000 ) second order semi-discrete central scheme in dimension by dimension 



formulation with a new two-dimensional approach introduced here and applied in numer- 
ical simulations for two-phase, two-dimensional flows in heterogeneous formations. This 



semi-discrete central scheme is based on the ideas of Rusanov's method (Rusanov (1961)) 



using a more precise information about the local speeds of wave propagation computed at 
each Riemann Problem in two-space dimensions. We find the KT dimension by dimen- 
sion has a much simpler mathematical description than the genuinely two-dimensional 
one with a little more numerical diffusion, particularly in the presence of viscous fingers. 



Unfortunately, as one can see in Abreu et al. (2007), the KT with the dimension by di- 



mension approach might produce incorrect boundary behavior in a typical geometry used 
in the study of porous media fiows: the quarter of a five spot. This problem has been 
corrected by the authors with the new semi-discrete scheme proposed here. We conclude 
with numerical examples of two-dimensional, two-phase fiow associated with two distinct 
fiooding problems: a two-dimensional fiow in a rectangular heterogeneous reservoir (called 
slab geometry) and a two-dimensional fiow in a 5-spot geometry homogeneous reservoir. 

1. INTRODUCTION 

We consider a model for two-phase fiow, immiscible and incompressible displacement in 
heterogeneous porous media. The highly nonlinear governing equations are of very prac- 



tical importance in petroleum engineering (Peaceman, 1977 Chavent and Jaffre, 1986) 



see also (Furtado and Pereira 2003) and the references therein for recent studies for the 



scale-up problem for such equations). 

The conventional theoretical description of two-phase fiow in a porous medium, in the 
limit of vanishing capillary pressure, is via Darcy's law coupled to the Buckley-Leverett 
equation. The two phases will be referred to as water and oil, and indicated by the 
subscripts w and o, respectively. We also assume that the two fiuid phases saturate the 
pores. With no sources or sinks, and neglecting the effects of gravity these equations are 



see Peaceman (Peaceman, 1977)): 



V ■ V = 0, 



-A(s)i^(x)Vp, 
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ds 
di 



+ V-(/(s)v) = 0. 



(2) 



Here, v is the total seepage velocity, s is the water saturation, K{x.) is the absolute 
permeability, and p is the pressure. The constant porosity has been scaled out by a 
change of the time variable. The constitutive functions A(s) and f{s) represent the total 
mobility and the fractional flow of water, respectively. 

We are concerned with numerical schemes for solving scalar hyperbolic conservation 
laws arising in the simulation of multiphase flows in multidimensional heterogeneous 
porous media. These schemes are non-oscillatory and enjoy the main advantage of 
Godunov-type central schemes: simplicity, i.e., they employ neither characteristic de- 
composition nor approximate Riemann solvers. This makes them universal methods that 
can be applied to a wide variety of physical problems, including hyperbolic systems of 
conservation laws. The two main classes of Godunov methods are upwind and central 
schemes. 



The Lax-Friedrichs (LxF) scheme (Lax, 1954) is the canonical first order central scheme, 
which is the forerunner of all central differencing schemes. It is a non-oscillatory scheme 
based on piecewise constant approximate solution and it also enjoys simplicity. Unfortu- 
nately the excessive numerical dissipation in the LxF recipe (of order 0{AX'^/At)) yields 
a poor resolution, which seems to have delayed the development of high resolution cen- 
tral schemes when compared with the earlier developments of the high resolution upwind 
methods. Only in 1990 a second order generalization to the LxF scheme was introduced 



in (Nessyahu and Tadmor, 1990). They used a staggered form of the LxF scheme and re- 



placed the first order piecewise constant solution with a van Leer's MUSCL-type piecewise 



linear second order approximation (Van Leer, 1979). The numerical dissip ation present 
in this new central scheme has an amplitude of order 0{AX'^ / At), (see (Abreu et al. 



2004c[ ), dAbreu et al.]|2004b[ ), ( |Abreu et al.]|2004aD , ( |Abreu et al.]|2006| ), for recent studies 



in three phase flows using the Nessyahu and Tadmor (NT) central scheme). In spite of 
the good resolution obtained by the Nessyahu and Tadmor scheme, much higher than in 
the first order LxF scheme, there are still some difficulties with small time steps which 
arise, e.g. in multiphase flows in highly heterogeneous petroleum reservoirs or aquifers. 
To overcome this difficulty, we can use a semi-discrete formulation coupled with an appro- 
priate ODE solver retaining simplicity and high resolution with lower numerical viscosity, 
proportional to the vanishing size of the time step At. Both LxF and NT schemes do not 



admit a semi-discrete form; see (Kurganov and Tadmor, 2000) for a detailed description of 



the one-dimensional Kurganov and Tadmor central scheme which is the first fully discrete 
Godunov- Type central scheme admitting a semi-discrete form. 



We compare the Kurganov- Tadmor (KT) two-dimensional (Kurganov and Tadmor 



2000) second order semi-discrete central scheme in dimension by dimension formulation 



with a genuinely two-dimensional approach applied in numerical simulations for two- 
phase, two-dimensional flows in heterogeneous formations. We flnd the KT dimension 
by dimension has a much simpler mathematical description than the genuinely two- 
dimensional one adding only a little more diffusion, particularly in the presence of viscous 
fingers. Unfortunately, the KT with the dimension by dimension approach might produce 
incorrect boundary behavior in a typical geometry used in the study of porous media flows: 
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the quarter of a five spot. These resuhs are presented in (Abreu et al., 2007). This prob- 



lem motivated the authors to develop a genuinely two-dimensional formulation which is 



then presented in section (2.2 ). Although a similar two-dimensional formulation was avail- 



able in a early work (Kurganov and Petrova, 2001), ours was developed independently to 



deal with two-phase flows, immiscible and incompressible displacement in heterogeneous 
porous media. It shares the same general ideas with the work of Kurganov-Petrovna but 
differs in many technical details. 

This paper is organized as follows. In Section [2] we introduce the model for two- 
phase flows, immiscible and incompressible displacement in heterogeneous porous media. 



In Section 12.21 we discuss the mathematical formulation for the KT central scheme in 
dimension by dimension approach and in a genuinely two-dimensional one. In Section [3] 
we will present some numerical results when we apply the KT central differencing scheme 
with both approaches mentioned above to porous media flows. 

2. NUMERICAL APPROXIMATION OF TWO-PHASE FLOW 

2.1. Operator splitting for two-phase flow. An operator splitting technique is em- 
ployed for the computational solution of the saturation equation ^ and the pressure 
equation ([T]) which are solved sequentially with distinct time steps. This method has 
proved to be computationally efficient in producing accurate numerical solutions for two- 
phase flow ( [Douglas et al. 1997). 

Typically, for computational efficiency, larger time steps are used to compute the pres- 
sure ([T|. The splitting allows time steps used in the solution of the pressure- velocity 
equation that are longer than those allowed in the convection calculation (|2]). Thus, we 
introduce two time steps: Ate used to compute the solution of the hyperbolic problem, 
and Atp used in the pressure- velocity calculation such that Atp > Ate- We remark that 
in practice, variable time steps are always useful, especially for the convection micro-steps 
subject dynamically to a CFL condition. 

For the global pressure solution (the pertinent elliptic equation), we use a (locally con- 
servative) hybridized mixed finite element discretization equivalent to cell-centered finite 
differences (Douglas et al. , 1997), which effectively treats the rapidly changing permeabili- 



ties that arise from stochastic geology and produces accurate velocity fields. The pressure 
and the Darcy velocity are approximated at times = mAtp, m = 0,1,2, ... . The alge- 
braic system resulting from the discretized equations can be solved by a preconditioned 
conjugate gradient procedure (PCG) or by a multi-grid procedure ((Douglas et al. , 1997[) ). 
We use high resolution numerical central scheme (see (Kurganov and Tadmor, 2000)) 



for solving the scalar hyperbolic conservation laws arising in the convection of the fluid 
phases in heterogeneous porous media for two-phase flows - we will discuss the application 
of these schemes for two-phase flows in Section |3] Theses methods can accurately resolve 
sharp fronts in the fluid saturations without introducing spurious oscillations or excessive 
numerical diffusion. 

The saturation equation is approximated at times = t"^ + nAt^ for t"^ < < t"^~^^ 
that take into account the advective transport of water. We will write to indicate the 
time step and t'^'^^ to indicate + Ate- 

A detailed description of the numerical method that we employ for the solution of Eqs. 
([T])-([2]) can be found in Douglas et al. (1997) and references therein (see also Abreu et al. 
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(2006, 2004c) and Abreu et al. (2008) for applications of the operator splitting technique 



for three phase flows that takes into account capillary pressure (diffusive effects)). 
i?emark: To simplify notation, we denote: 

• NTld for one-dimensional NT scheme; 

• NT2d for two-dimensional NT scheme; 

• KTld for one-dimensional KT scheme; 

• KTdxd for the KT scheme with dimension by dimension approach and 

• SD2-2D for our two-dimensional approach. 



2.2. TWO SPATIAL DIMENSIONS SECOND ORDER SEMI-DISCRETE CEN- 
TRAL SCHEME. In this section, we will develop a two-spatial dimension second order 



semi-discrete central scheme (SD2-2D) based on the ideas of 


Lax 


(1954 


); 


Rusanov 


1961) 


Nessyahu and Tadmor 


1990 ) ; Kurganov and Tadmor 


(2000 


) and . 


fiang and Tadmor 


(1998 



which are then applied in numerical approximation of the scalar hyperbolic conservation 
law modeling the convective transport of the fluid phases in two-phase flows and its cou- 



pling with lowest order Raviart- Thomas (Raviart and Thomas, 1977) locally conservative 



mixed finite elements for the associated elliptic (velocity-pressure part) problem (See 
Raviart and Thomas (1977)). We summarize below the basic ideas of the construction of 
SD2— 2D numerical scheme: 

• The Lax-Friedrichs method in two-spatial dimensions LxF2D written in the REA 



algorithm setup (See Jiang and Tadmor (1998)) will be used to obtain the two 



dimensional Rusanov's method SD1-2D. We follow the same procedures presented 
in 



Kurganov and Tadmor (2000) in one spatial dimension. 
The new SD2-2D numerical scheme will then be obtained as a second order gen- 
eralization of the SD1-2D. This second order precision is achieved approximating 
the solution with piecewise linear functions. 



2.3. The staggered non-uniform grid of the SD2-2D central scheme. We begin 
then extending the LxF2D to obtain the SD1-2D following the same procedures for one 
dimensional problems. First, we define the non-staggered and the staggered grids of 
retangular cells used in the LxF2D. The points {xj,yk,t'^) are defined as follows. 

Xj = j-Ax, j = ...,-1,0, 1,... 
yk = k-Ay, k = ..., -1,0,1, .. . 

We denote the cells of the non-staggered grid by Ij^k '■= (a;j-i/2, a^i+1/2) x (yfc-1/2, l/fc+i/2)- 
Its area Aj^k = Ax- Ay = (xj+1/2— Xj_i/2)-(?/fc+i/2— l/fc-i/2)- The time step of the convective 
equation ^ is Ate = ^^^^ ~ t"^- The staggered grid is obtained moving the cells Ax/ 2 
to the right and Ay/2 upward These staggered cells will be denoted by Ij+i/2,k+i/2 '■= 
{xj,Xj+i) X {yk,yk+i)- Its center is the point (xj+1/2, yfc+1/2), where Xj+1/2 = xj + Ax/2 
and yk+1/2 = yj + Ay/2. The Figure [l] illustrates the non-staggered and the staggered 
grids showing the cells Ij^k e Jj+3/2,fc+i/2 as examples of non-staggered and staggered cells, 
respectively. 
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Figure 1. The LxF2D uniform grid. The cells 1^ ^ of the original non- 
staggered grid are limited by the solid lines and the cells /j+3/2,fc+i/2 of the 
staggered grid are limited by the dashed lines. 



The scalar hyperbolic conservation law ^ can be written as 



(3) 



where ^ = ^{x,y,t) and % = 'hj{x,y,t) denote the x and y components of the velocity 
field V. The cell averages at time are 



AXAY 



J+2 / (-H-j 



J 2 2 



s{x, y, f^) dxdy. 



(4) 



The solution s{x,y,t'^) of the probl em ([2| at time is approximated using piecewise- 
linear MUSCL-type interpolants (See Van Leer (1979)). 



SU^, y) = si, + is,)lk ■ {x - xj) + iSy)l, ■ {y - y,), 



(5) 



where a;j_i/2 < x < e yk-1/2 < 2/ < Uk+i/i- The second-order resolution is guaran- 

teed if the so-called vector of numerical derivative {Sx)j^k ('^y)j,A: satisfy 



{SyYj^k 



ds 
dx 

ds 
dy 



+ 0(AX); 



x=Xj,y=yk,t=i"- 



+ 0{AY), 



(6a) 
(6b) 



x=Xj,y=yk,t=t>^ 
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These numerical derivatives are computed using the MinMod hmiter 

\ Ax 2Ax Ax J 

iSy)ik = MM^i-{:s;,_„:s;„:s;,,j 

V Ay 2Ay ' Ay J ^ ^ 

Here A denotes the centered difference, ASj^-^^^^ k ~ ^'j+i,k ~ ^'j,k ^^'^ ^^e paramter 6 G 
[1, 1.8] has been chosen in the optimal way in every numerical example with 6 = 1.8 
beeing the less dissipative limiter. The minmod limiter ([T]) guarantees the nonoscillatory 
property and the second-order accuracy. The reconstruction given by Equations (|5])-([7]) 



also retains conversation, i.e.. 



Slk{x,y)dxdy = S]u. (8) 



Remark: We notice that if {Sx)j f, and {Sy)ji^ are equal to zero, then we will get the 
first-order two-dimensional semi-discrete scheme SD1-2D. Otherwise, we will obtain the 
second-order two spatial dimensions semi-discrete central scheme SD2-2D. 

We consider the model of hyperbolic conservation laws given by Equation ^ with cell 
averages as in Q and the two-dimensional piecewise linear reconstruction defined in ([s]) 
and ^ with the conservative property ([s]). Our goal is to compute an approximated 
solution Sj^kiim + Ate) in the original grid at the next time step. To this end, we apply 
the Godunov's algorithm REA. To solve this problem, we integrate the conservation law 
over some control volumes that we need to specify. 

Constructing the staggered nonuniform grid: Kurganov and Tadmor developed 



the KTID scheme along the lines of NTID (See (Kurganov and Tadmor 2000)). The 
nonuniform staggered grid in the KTID was constructed directly from the staggered 
uniform grid of NTID with additional information on the local speeds of wave propagation. 
In a similar way the nonuniform staggered grid in the SD2-2D scheme is defined from 
the uniform staggered grid of the NT2D scheme as follows. 

(1) We set Cj^k = [xj-i/2,Xj+i/2] x [?//,•- 1/2, l/fc+1/2] to denote the cells of the original 
non-staggered grid; Cj+i/2^k+i/2 = [^^j^^jj+i] x [yk,yk+i] to denote the cells of the uniform 
staggered grid. We start with a piecewise constant approximated solution Sjj^ over the 
original cells Cj^k- 

(2) Next, we move to the staggered uniform grid with cells given by 6*^4.1/2,^+1/2 • 

(3) The NT2D scheme computes four averaged solutions at time step + At^ (See 
the hachured regions in Figure [2] corresponding to the cells Cj_ 1/2,^-1/2, C'j+i/2,fc-i/2, 
Cj-i/2,fc+i/2 and Cj+i/2,fc+i/2 on the staggered grid). These averaged staggered solution 
are then properly projected back onto the original non-staggered grid to obtain the desired 



solution (See Jiang and Tadmor (1998)). 
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Repeating the same ideias presented by Rusanov in his modification of Lax-Friedrichs' 
method, we compute the local speed of propragation at each Riemann Problem. These 
local speeds define new non-uniform cells where the evolution step will take place. 

Computing the local speed of propagation: we begin with the cell Cj i/2,fc 1/2 
to find the local speeds at the following Riemann Problems: 

(1) Y direction: 

(a) 

5^_i,fc-i, a;j_3/2 < a; < x^.i/a, yk-z/2 <y < 
^j-i,k^ ^j-^l-i < a; < a;j_i/2, yk~i/2 <y < yk+i/2- 
The local speed aj_^ fc-1/2 defines the following points: 

Pi = (xj^i, yk^i/2 - ay_-^ ^,_-^^^Atc), 



sketched in Figure|2| We denote the distance between them by A?/j_i fc_i/2 

(b) 



^j-if-i <x < Xj+i/2, yk-3/2 <y < yk-1/2 
^j,k^ <x < Xj+i/2, yk-1/2 <y < yk+i/2- 

The local speed aJfc_i/2 defines the points 

P3 = {Xj, yk-1/2 - aj;^_i/2^tc), 

Pi = {xj, yk-1/2 + aJ^fc_i/2Atc) 

shown in Figure |2jand the distance between them is Ayj^k_i/2 = 2a!'. ^_^i^Atc- 
(2) X direction: 



(a) 



S^-i^fe-i, a;j-3/2 < a; < Xj_i/2, l/fc-3/2 <y < yk-1/2 



Sj,k-i^ Xj-1/2 <x< Xj+i/2, yk-3/2 <y < yk-1/2- 
The local speed ci^_i/2 k-i defines the points 

P5 = {xj-1/2 - a^_i/2,k-Atc, yk-i), 
Pe = {xj^i/2 + aJ„i/2,fe_iAtc, yk-i) 

also shown in j2j and Axj^i/2^k-i = '^(^j_i/2 k-i^^c is the distance between 
them. 

(b) 

^l-i,k^ Xj-3/2 <x < Xj^i/2, yk-1/2 <y < yk+1/2 
^j,k^ Xj-1/2 <x< Xj+i/2, yk-1/2 <y < yk+1/2 ■ 
The local speed aJ_i/2A: defines the points 

P7 = {xj-1/2 - aJ„i/2,fcAtc, yk), 
Ps = {xj-1/2 + a^j_i/2,k^tc,yk) 



F. Furtado, F. Pereira and S. Ribeiro 



shown in Figure [2j and AXj_i/2,fc = fc^^c denotes the distance between 

them. 



Given these four local speed of wave propagation a^ -^ ^ 



we can define the Region I (also represented by Rj-i/2,k~i/2) as follows: 
> Region I: 



y X X 

/2' '^j,A;-l/2' ^j-l/2,k-l ^ '^j-l/2,fc' 



-Rj-i/2,fc-i/2 := 



[Xj-l/2 - &J_i/2,fc-l/2^^c, + ^J-l/2,fc-l/2'^y ^ 

[l/fc-1/2 - ^|-l/2,fc-l/2^^c, 2/fe-l/2 + ^J-l/2,/c-l/2^^c] 



where feJ_i/2,fc-i/2 — max{aj_i/2,fc, aj_i/2,fc-i} e 

^i-l/2,A:-l/2 "^^^{'^j,A:-l/2' ^j-l,fc-l/2}- 



Figure |3| shows the new cell Rj-i/2,k-i/2 of the new staggered non-uniform grid. 
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Figure 2. SD2-2D: the construction of the two-dimensional grid 



We repeat the same procedures with the cell Cj_ 1/2,^+1/2 to define the Region III in 
terms of the local speed of propagation aJ_]^^^]^/2' '^jfc+1/2' '^j-i/2fc+i ® ^j-i/2fc- These 
local speeds determine analogously the points q\ a q% sketched in Figure [3] and define the 
new Region III: 
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> Region III: 

Rj-l/2,k+l/2 '■- 



[Xj^l/2 — ^j_l/2,A:+l/2^^c, Xj-1/2 + ^j-l/2,fc+l/2^^c 
[yk+1/2 - &y_l/2,fc+l/2^^c, yk+1/2 + h)_i 



i-l/2,/c+l/2 



At, 



where &J_i/2,fc+i/2 := max{aj_i/2,fc, aj-i/2,fe+i} and 

^j-l/2,fc+l/2 '"^^'^{^\k+l/2^ ^^j-l,k+l/2} ■ 



k+1/2 



' k-1/2 





i q 
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Figure 3. Regions I and III 

Following tliese same procedures witli tlie staggered cells Cj+i/2,k-i/2 and Cj+i/2,fc+i/2 we 
define Regions VII and IX shown in Figure |2.3[ We will denote by Group A the set of 



Regions I, III, VII and IX. Finally, to finish the construction of the new non-uniform 
staggered grid, we need to define: 

• Four cells which lie in the empty spaces between the regions of Group A. This new 
set will be denoted by Group B. 

• A central region where the solution is smooth. 

The definitions of the cells of Group B will determine the central region. This central 
region will not be a rectangle, but a set of retangles. There are many ways to define the 
cells of Group B. Our definition will be as follows (See Figurajs]): 
> Region II: 

Rj-l/2,k '■= [Xj-1/2 - cJ_i/2,fcAtc, Xj_i/2 + C^_i/2,k^^c] X 

[yk-1/2 + b^j_^/2,k-l/2^^c, yk+1/2 - ^J-l/2,fc+l/2^^c] 



where 



Cj-1/2,A: ^ax{6J_]^/2,fc-l/2) ^J-l/2,fc+l/2} 
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Vk+I 



y k+1/2 



yk 



y k-1/2 



^k-l 



^j-1 ^ j-1/2 ''j ''j+1/2 ^j+1 

Figure 4. Regions I, III, VII and IX 

The Regions VI and VIII can be obtained analogously As soon as the cells of Group A 
and B are determined, the central region is automatically defined: 

> Region V: 

Rj,k '■= [^j-l/2 + cj_i/2,A:^^c, - C^^-^/2,k^'l^c] X 

[yk-1/2 + dl^^y^Atc, yk+1/2 - rfj^fc+i/2^tc]. 

We also would like to emphasize that our choice for these regions does not introduce 
more numerical diffusion. We will call BR the black rectangle that can be seen in Figure 
[sjand we notice that, by construction, its area is proportional to (Atc)^- 

2.4. The new SD2-2D central scheme using Algorithm REA. After defining the 
new control volumes performed in the section above, we are now able to develop our new 
SD2-2D central scheme following the REA algorithm. 

Reconstruction step: We suppose that we know an approximated solution constant 
in each cell at time step as in Equation Q. This approximated solution is then 
reconstructed as a piecewise bilinear polinomial as defined in Equations ^ and (|6]). 
Evolution step: Let D represents one of the nine regions defined above Rj±i/2,k±i/2, 
Rj,k±i/2, Rj±i/2,k, the central region Rj^ or the black rectangle BR. We will denote by 
the part of region D inside the non-staggered cell Ij^k and by , the part of region 
D outside the cell fc. 

We integrate the conservation law ^ in the control volumes D x [t^, t^ + Atc] to obtain 
an approximate averaged solution w'^~^^{D) at the next time step, in each cell D of the 
staggered non-uniform grid. 
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' k+1 



' k+1/2 
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N-1 



^ j-1/2 



Figure 5. SD2-2D: The two dimensional semi-discrete central scheme 
SD2-2D: construction of the non- uniform staggered grid. 

Projection step: These averaged solutions w'^^^{D) are then reconstructed as piece- 
wise bilinear polynomials w'^~^^{x,y) in each of the ten regions D. These new reconstruc- 
tions are then projected back onto the original grid of uniform non-staggered cells, 

— ^ / w^+\x,y)dxdy. (9) 
Ax Ay 



S 



■K+l 



The new reconstructions w'^^^{x,y) are defined analogously as in Equation (|5]). For 
instance, for Region Rj+i/2,k, 

{x,y) e Rj+i/2,k- 

The numerical derivatives ('W^a;)j+i/2 ^"^^ ('"^y)j+i/2 fc satisfy the conditions 





n _ 


dw 




-1/2, fc 


dx 


K- 


n _ 


dw 




-1/2, fc 


dy 



+ 0{Ax) 



(11) 

(12) 



+ 0{Ay)- 

in order to guarantee the second order approximation. Also, the reconstruction Wj^y2 fc(^' ^) 
retains the conservation property ([s]). We remark that this is a theoretical step and it 
will not be necessary to compute these numerical derivatives in the final semi-discrete 
formulation. 
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This completes the construction of our totally discrete central scheme in a rectangular 
grid. It is very laborious to write a totally discrete version of this central scheme. Instead, 
we will proceed directly to our semi-discrete formulation. In order to to this, we compute 
the following limit when At^ 0, 



lim + ^^^^ " ^^^^^^^ - 



Atr 



w 



lim - — ■ - — — 
At^^o Ate Ax Au I — / R 




p,fc+l/2 



:_Li /n ^ 7, 1 /n -IT /o ^ R„ u 



p=j±l/2 P,fc-l/2 



p=i±l/2 ^V,*; 



+ Y w]^\x,y)dxdy+ w']^\x,y) dx dy 

q=k±l/2 -^^ti "^^J'fe 

+ / w^j+i/2,fc_i/2(a;>l/)'^a;rf?/ - (Ax A?/)S'j-fe(t) 



(13) 



The conservation property of the reconstructions u)j^^ in he regions D results 



w 



■K + l, 



X, y) dx dy = \D\ ■ w'^^'^{D) 



(14) 



D 



where w'^^^{D) is the averaged solution in region D. Note that, by reconstruction, the 
area of regions I, III, VII and IX are proportional to (Atc)^, that is. 



Ref 


^ion 


I: 


-Ri-l/2,fc-l/2 


= O 




Ref 


^ion 


III: 


-Rj-l/2,fc+l/2 


= O 




Ref 


^ion 


VII: 


-Rj+l/2,fc-l/2 


= O 


(Ate)^ 


Rej 


^ion 


IX: 


-Ri+l/2,fe+l/2 


= O 




Rej 


^ion 


RP: 


\RP\ = 0[{At,f) 





For example, considering Region IX, we conclude 



i+l/2,fe+l/ 



i/2i^^y) dxdy = 0{Atcf 



7 + 1/2, fc+l/2 



lim ^ 

Atc^O At, 



w 



c JR 



K + 1 

i+l/2,A:+l/2 



[x, y) dx dy = 0; 



(15) 



-j + l/2,fc+l/2 



Our goal is to obtain the semi-discrete formulation of this new central cheme. Therefore, 
the Equation ( 15 ) show that we do not need to computed the averaged soltions over the 
Regions I, III, VII and IX. Note that, by simetry, we only need to compute the solutions 
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over the Regions VI, VIII, V and the Black Rectangle. For the Region Rj+i/2,k, we obtain: 



j + l/2,fc 



+ 1 I p + 

"^i+l/2,fc ■ l-^j+l/2,fc 



dxdy 



+ 0((Ate) 



(16) 



Analogously, we compute the averaged solution over Region R^f._^_i/2- 



(17) 



Note that the solution has no discontinuities inside Rj^k- So, it isn't necessary to recon- 
struction as a piecewise bilinear polynomials. The averaged solution is 



/ w'J^\x, y)dxdy = w]^^ ■ \Rj^k\ 



;i8) 



Substituting the Equations (15), (16), (17) and (18) in Equation (13), we obtain: 



d_ 
It 



Sj,kit) 



'i-l/2,A:- 



C 



j+l/2,k 



Ax '^.~^Mt + ^Q+ 



Wj+l/2,k{t + Ate) 



lim 

d^ d^ 

+ ^4^%fc_l/2(t + At,) + ^^^w,^k+l/2{t + At, 



Ay 



^],k+l/2 + ^\k-l/2 ^j+l/2,k + ^j-l/2,k 



Ay 



Ax 



Wj^kit + Ate 



+ ( -^""l^kit + Ate) - ^S,^k{t + Ate 



At, 



(19) 



For the final formulation, we have to compute the averaged solution over the non- 
uniform staggered grid. 



^j+l/2,k-> ^j,k+l/2 ^ ^j,k ■> 



To this end, we integrate the conservation law (??) over the control volumes 



Rj+l/2,k X [CjC + Ate], Rj,k+l/2 X [C>^m + At,] C Rj^k X [t^,t^ + At, 
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respectively. Therefore, 



w 



K+l 



1 



\R 



1 



'j+l/2,k\ JRj+i/2,k 



s{x, y, dx dy 
S'^{x, y) dx dy 



1 f [ \^( 



''v{x,y,T) ■ f{s{x,y,T)) 



(20) 



+ 



d_ 
dy 



'^v{x,y,T) ■ f{s{x,y,T)) 



dx dydr. 



Let us denote the double integral by Intg and the flux integral by Int/. The integral Int^ 
is computed analytically. 



1 



+ 



4 
At 



(21) 



^j+i/2,it-i/2 ^j+i/2,fe+i/2 ■ [i^y)j,k i^y)j+i,k] ■ 



To compute the flux integral, Int/, we first denote the limits of region i?j_|_i/2,jt as follows: 



Atr 



Xj+l/2 — Cj^i/2,k 
Xj+1/2 + Cj+i/2,fcAtc 
l/fc+1/2 + ^J+l/2,fe-l/2^*c 
Z/fe+1/2 - ^'J+l/2,fe+l/2^^c 



Using the Calculus Fundamental Theorem together with the trapezoid rule, we obtain 
Int/ 



2Ax 



—f 



h{b, d, t) f{s{b, d, r)) - 't;(a, d, r) /(s(a, d, r)) 

dr 

yv(b, d, r) f(s(b, d,r))- Mb, c, r) f(s(b, c, r)) 

dr (22) 



+''v{b, c, r) /(s(6, c, r)) - ^^(a, c, r) f{s{a, c, r)) 

2Ayj+i/2,fc ,A« 
+%(a, r) /(s(a, d, r)) - %(a, c, r) /(s(a, c, r)) 



If the CFL condition 



/At, At, A 1 

max I — max I v / (s)|, — max |% / (s)| 1 < - 



(23) 



holds and since the functions ^ f{s{x, y, r)) and % f{s{x, y, r)) are computed away from 
the discontinuities then they are differential functions of r and therefore, the time integral 
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can be approximated using the middle point rule. Denoting t'^+^Z^ ;= t + Atc/2, we obtain: 



Int; 



^ Xb, d, f{S{b, d, - ""via, d, t-+i/2)) 



dr 



ay 

H • 

1/2 ^j+l/2,fc-l/2) 

+^t;(a, d, t-+i/2)) _ y^^a^ t'+^Z') /(S(a, c, rfr (24) 

where ax = At J AX and ay = At J AY. 

The midpoint values are computed using the Taylor expansions and the conservation 
law For instance, 

' ^(a,d,r+i/2) := S{a,d,t) 

''Xa, d, t) f{S{a, d, t))] (Xa, d, t) f{S{a, d, t)) 



S{a,d,t) := S^^^^k - Ax (5^.)i+i,fc(^ - "xcJ+i/2,fc)-Ay {Sy)]^^,^(^^ - tty&J+i/2,fc+i/2)- 
As the time step Ate goes to zero, the limit 

^I^^Q^i^^^^^'^^^^'^) ~ ^j+i,k Y~ ('^^''i+i.fc 2~ ^^y^J+'^^>' '~ '^j+i/2,fc-i/2- (25) 
These are called the intermediate values and their general form is 

'^i+l/2,fc+l/2 ^ ^ j+l/2±l/2,k+l/2±l/2 i ('S'z)j+l/2±l/2,fc+l/2±l/2(^i+l/2 " 2^j+l/2±l/2) 

i-;^- ('S'y)j+l/2±l/2,fc+l/2±l/2(l/fc+l/2 " l/A;+l/2±l/2) (26) 



We notice that the cell averages w'^\^ are obtained analogously to (20). 

And also, w'r^lj^ ,^ = w;^l^^_,^, e w;J^\^^ = <^+i/2_r 

Substituting all these cell averages in time step + Ate into Equation ( 19 ) and com- 
puting the limit when At^— >0, we obtain the second order central scheme in semi-discrete 
formulation: 

d^ ... -^J+l/2,fc(^) --^J-l/2,fc(^) -^Ifc+l/2(^) --^Ifc-l/2(^) 

dt^^'^^> AX AT ' 
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where the numerical fluxes are 



Hj+l/2,kit) -^yVj+i/2,k+l/2{t) /(5'+^i/2,fc+l/2(^)) + /('S'j+i/2,fc+l/2(^)) 

'^j+l/2,A:('^) " ^j+l/2,ki'^) 



Cj+l/2,fc 



^^+1/2^ = 4|Vl/2,fc+l/2W [/('5,-++i/2,fc+i/2W) + /('5^+i/2,fe+l/2W) 
+ ^-l/2,.+l/2(t) [/(5ri/2,.+l/2W) + /(^;-:i/2,.+l/2W)' 



"i,fc+l/2 



s 



j,k+l/2 



(t) - s. 



j,k+l/2 



it) 



(28a) 



(28b) 



If the numerical derivatives are equal to zero then we obtain the two-dimensional Ru- 
sanov's central scheme in semi-discrete formulation. 



RusaJ^_i/2,fc(^) 



1 



"^^^+1/2,^+1/2 (t) /(5,-+l,,(t)) + /(^,,fe(t)) 



+ Vl/2,fc-l/2(t) [/(5', + l,fc(t)) + fiSj,kit)) 



Cj+l/2,fc 



R'UsaJfc^^/2W =7r^j+i/2,fc+i/2(t) fiSj,k+iit)) + fiSj^kit)) 



+ ^,-_l/2,fc+l/2(t) [/(5,,fc+l(t)) + f{S,4t)) 

Sj,k+i{t) — Sj^kit) 



i,fc+i/2 



(29a) 



(29b) 



This new two-dimensional semi-discrete central scheme with the numerical fluxes given 



by (28) or (29) comprises a system of ordinary differential equations. To solve this system, 



we use the explicit second order Runge-Kutta method. 



2.5. The velocity field. Finally, to complete the description of the genuinely two- 
dimensional KT scheme, we have to define the velocity field. The velocity is defined 
at the vertices of the cells. We can not use directly the velocity field from the Raviart- 
Thomas space as we did in the dimension by dimension approach. Instead we will use a 
bilinear interpolation of it preserving the null divergence necessary for the incompressible 
flows. For instance, to compute the value of "^^+1/2,^+1/2 on the vertex (xj+1/2, i/fc+1/2) at 
some time step t™, we have to use all the four cells which share this vertex. 



SD2-2D 
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f^j+l/2,fc+l/2 



2 

1 rl 



+ 



2(2' 



^rj+l/2,fe ~ 'Vlj-l/2,k) 
+ 7ji'"rj+l/2,k+l — 'Vlj-l/2,k+l) + (:^('^Vj+3/2,A: ~ ^j+l/2,fc) + 

'^rjr-+l/2,A; ~ ^j-l/2,fc + "^^r j+l/2,fc+l ~ ''^j-l/2,A;+l 
+'^r-j+3/2,A; ~ ^j+l/2,fc + ''^r j+3/2,fc+l ~ ^lj+l/2,k+l^ 

3. TWO-DIMENSIONAL NUMERICAL EXPERIMENTS 

We present and compare the results for numerical simulations of two-dimensional, two- 
phase flow associated with two distinct flooding problems using the KT scheme with the 
dimension by dimension approach (KT dxd) and the genuinely two-dimensional formu- 
lation (SD2-2D). The first problem is a two-dimensional flow in a rectangular heteroge- 
neous reservoir (called slab geometry) with 256 m x 64 m in size, and the second is a 
two-dimensional flow in a 5-spot geometry homogeneous reservoir having 64 m x 64 m. 

In the 5-spot geometry homogeneous reservoir simulation we used two distinct uniform 
five-spot well configurations intended to illustrate different flow patterns, with parallel 
and diagonal grid orientations, and boundary behavior. 

In all simulations the reservoir contains initially 79% of oil and 21% of water. Water is 
injected at a constant rate of 0.2 pore volumes every year. 

The fractional volumetric flow, the total mobility, and the relative permeabilities are 
assumed to be: 

k 

X{s) 



fis) 



and 



Xis) 



Sj.Q 



1^0 



-^(s-s ? 



where Sro = 0.15 and Srw = 0.2 are the residual oil and water saturations, respectively. 
The viscosity of oil and water used in all simulations are fio = 10.0 cP and fio = 0.05 cP. 
For the heterogeneous reservoir studies we consider a scalar absolute permeability field 



K{x.) taken to be log-normal (a fractal field, see (Glimm et al 





1993 


) and ( 


Furtado 



and Pereira, 2003) for more details) with moderately large heterogeneity strength. The 
spatially variable permeability field is defined on a 256 x 64 grid with three different 
values of the coefficient of variation C„ (standard deviation) /mean: 0.5, 1.0, and 2.4. 

The boundary conditions, injection and production specifications for two-phase flow 
equations ([T])-([2])) are as follows. For the horizontal slab geometry, injection is made 
uniformly along the left edge of the reservoir and the (total) production rate is taken to 
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be uniform along the right edge; no flow is allowed along the edges appearing at the top 
and bottom of the reservoir in the graphics (Figures [6| [7| and [8|. 

In the case of a five-spot flood with diagonal grid (Figure joj), injection takes place 
at one corner and production at the diametrically opposite corner; no flow is allowed 
across the entirety of the boundary. In the case of a five-spot flood with the parallel grid, 
injection takes place at two corners (diametrically opposite), say left down and right up, 
and production in the diametrically 'off corners', say right down and left up. 



SD2-2D 
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Figure 6. Water saturation surface plots for two-phase flow in a two- 
dimensional heterogeneous reservoir having 256 m x 64 m, with the coef- 
ficient of variation = 0.5 and viscous ratio 20. From top to bottom:!) 
(KT dxd) scheme with 256 x 64 grid; 2) (KT dxd) scheme with 512 x 128 
grid; 3) (KT two) scheme with 256 x 64 grid. 
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Figure 7. Water saturation surface plots for two-phase flow in a two- 
dimensional heterogeneous reservoir having 256 m x 64 m, with the coef- 
ficient of variation CV = 1.2 and viscous ratio 20. From top to bottom:!) 
(KT dxd) scheme with 256 x 64 grid; 2) (KT dxd) scheme with 512 x 128 
grid; 3) (KT two) scheme with 256 x 64 grid. 



SD2-2D 
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Figure 8. Water saturation surface plots for two-phase flow in a two- 
dimensional heterogeneous reservoir having 256 m x 64 m, with the coef- 
ficient of variation CV = 2.2 and viscous ratio 20. From top to bottom:!) 
(KT dxd) scheme with 256 x 64 grid; 2) (KT dxd) scheme with 512 x 128 
grid; 3) (KT two) scheme with 256 x 64 grid. 
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16 32 48 64 16 32 48 64 



(c) NT2D: 128x128 clulas (d) KTdxd: 128x128 clulas 

Figure 9. Water saturation level curves for two-phase flow in a five-spot 
well configuration. The SD2-2D scheme was used in pictures in the left 
column and the KTdxd was used in pictures in the right column. 



SD2-2D 
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3.1. Analyzing the Numerical Results. Conclusions. The Figures |6} [7] and [8] refer to 
a comparative study for the KT dimension by dimension and a genuinely two-dimensional 
KT schemes showing the water saturation surface plots after 350 days of simulation for 
three different values for the strength of the heterogeneity of the fractal permeability field, 
CV = 0.5,1.2 and 2.2. 

Note that the genuinely two-dimensional KT scheme gives a more accurate solution 
than the solutions computed by the KT dimension by dimension scheme for the same grid. 
In fact we observe that the KT dxd is only comparable in accuracy with one degree of 
refinement (see Figures [6| [T] and |8]). The better accuracy of the genuinely two-dimensional 
approach is due to a more precise computation of the genuinely two-dimensional numerical 
fluxes, with respect to the one dimensional numerical fluxes in the dimension by dimension 
approach. 

In the case of a five-spot geometry homogeneous reservoir. Figure [9] (diagonal grid) 
shows the saturation level curves after 260 days of simulation obtained with KTdxd and 
SD2-2D schemes for two levels of spatial discretization. In this figure [9} the pictures on 
the left column are the results obtained with the SD2-2D scheme and the ones on the 
right were computed with the KTdxd scheme. In these Figures, the grid become finer 
from top to bottom, having 64 x 64 and 128 x 128 cells, respectively. 

It is clear that the KTdxd scheme (right column pictures in Figures [9] is producing 
incorrect boundary behavior. Moreover as the computational grid is refined (right column 
and bottom picture in Figure [9]) this problem seems to be emphasized indicating that the 
solution is not convergent. 

The KTdxd scheme uses numerical fluxes in the x and y directions which can be viewed 
as generalizations of one-dimensional numerical fluxes. We state that this type of approx- 
imation for the fluxes leads to the incorrect boundary behavior discussed above. This in- 
correct boundary behavior led us to develop a new genuinely two-dimensional KT scheme. 
The results obtained with this new scheme can be seen in the left column of Figure |9| It 
is clear that we have corrected the boundary behavior just changing the approach from 
dimension by dimension to our two-dimensional approach. This fact indicates (but does 
not prove) our idea that computing two-dimensional numerical fluxes using straight gener- 
alizations of one-dimensional numerical fluxes may produce incorrect numerical solutions. 
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